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CHAPTER I 


INTRODUCTION 

This research involves the use of deconvolution methods on gas chromato- 
graphic data to obtain an accurate determination of the relative amounts of 
each material present by separating the merged peaks. 

The underlying assumption implied throughout is that the area under the 
peaks is an accurate indicator of the amount of material present.^ Due to 
uncontrollable errors at every stage of operation, this is never exactly 
true. However, it is a consideration that must always be borne in mind in 
designing, building, and operating the gas chromatograph, its detector, and 
other related equipment. As peak areas are less dependent on operating 

conditions than other methods, they are preferred by most researchers for 

2 

accurate determinations. 

The promising new concept of function continuation in the frequency 
domain was introduced in striving to bring the art of deconvolution to its 
theoretical limit of accuracy. This proved only partially successful. 

The methods used here will have general applicability to any deconvolu- 
tion where the desired result consists of a number of well separated distrib- 


utions. 



Gas Chromatography 


Gas chromatography is an effective and widely used method for the 
separation of gases and volatile liquids or solids in the gaseous state. 

A small sample of the material being examined is injected into a stream 
of an inert gas such as nitrogen, hydrogen, carbon dioxide, argon, or helium 
which carries it into a column containing a suitable medium capable of retard- 
ing the flow, by varying degrees, of the individual components of the sample 
as they flow through the column. Differences in the time various components 
remain in adsorption or partition on the material in the column is again the 
factor which makes separation possible. The separated components then emerge 
from the column at discrete intervals (characteristic of each component) and 
pass through some form of detector. As a general rule, gas analyses are 
carried out on adsorption columns (gas/solid chromatography), while liquids 
and volatile solids are separated on partition columns (gas/liquid chromato- 
graphy).^ 

The data used in the present study were obtained from a gas/liquid 
chromatograph. The carrier gas was helium and the detector was a Varian 1400 
flame ionization detector. The sample compounds consisted of five xylenes. 

The flame ionization detector is a differentiating type, making the 
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area under the peak a suitable measure of the amount of material present. 

The flame ionization detector operates on the principle that the 
electrical conductivity of a gas is directly proportional to the concentration 
of charged particles within the gas. Effluent gas from the column is mixed 
with hydrogen and burned in air. Ions and electrons formed in the flame 



enter the electrode gap, decreasing the gap resistance, thus permitting a 
current to flow in the external circuit.^ 

The flame ionization detector is extraordinarily insensitive to air and 
water, making it especially suitable for the analysis of air pollutants or 
aqueous samples such as beverages, biological materials, and other liquids. 

For other materials the areas obtained must be multiplied by the proper 
correction constant to obtain the true proportions.^ 

The flame ionization detector has the widest linear range of any detector 
in common use.^ 

The analog output of the gas chromatograph was digitized by bringing its 

g 

output into the HP2100 minicomputer. Data were acquired at a 10 hz rate. 

Electronic digital integration was accomplished by accumulating the 

data points when an increase in signal was detected and terminating when the 

9 

signal again increased, which, incidentally, began a new count. 

Temperature programming was used on the substances making up this data. 
The temperature was increased at a linear rate different for each run.^^ 

The Data 

Data are given for five runs. The same sample composition, consisting 
of five test substances (see Figures 1 and 2), was injected each time. The 
graph of Run 3 was not included as nothing new of significance was illustrated 
Each run was taken at a different temperature. See Table 1 for specifications 
The interval between each data point was quite small compared to the 
standard deviation of the narrowest peak, thus making the data suitable for 
application of the fast Fourier transform. 
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Fig. 1. Data of the first two runs; (a) is data of Run 2, (b) is data of 


Run 1. 








TABLE 1 


GAS CHROMATOGRAPH SPECIFICATIONS 


Temperature Programming 


Run 1 

Run 2 

Run 3 

Run 4 

Run 5 

Beginning 

Temperature 

50° 

75° 

100° 

135° 

tn 

O 

O 

Rate of 
Change 

4°/mi n . 

A^^/min. 

4®/min. 

6°/min. 

6°/mi n . 


GX. Parameters 
10 '^ 

He 20 cc/mip. 

H 30 cc/min. 

2 

Air 300 cc/min. 
Injection 0.75y£ 22A 


Integrator 

TOO mv input = 1 v output 

Slope sens. 0.01 up; 1.0 down 

Zero delay 

Noise suppression 3 

Area threshold 1 

Shoulder control on (front) 
rear 0.1 



An inspection of the data points and its graph reveals that it is 
relatively smooth and free from noise. 

The tic marks, quite noticeable on the graph, were imposed on the 
data whenever the integrator encountered any increase in the signal. Removing 
the tic marks was quite easy. Graphs were drawn which included a few points 
on both sides of the tic marks. A smooth curve was drawn through these 
points. Intermediate points were then taken from the graph and replaced 
those of the tic marks on the data. The Interval was small and the data was 
digitized so that this method was sufficiently accurate for the purpose. - 
The baseline drift is negligibly small and the baseline will be taken 
at zero. 
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CHAPTER II 


FAST FOURIER TRANSFORM DECONVOLUTION 


The fast Fourier transform is a computational algorithm for calculating 
the discrete form of the complex Fourier series. (See Appendix B for a 
discussion of the discrete Fourier transform.) The addition theorem, shift 
theorem, convolution theorem, and other theorems familiar to users of the 
continuous Fourier transform with infinite limits may also be developed for 
the complex Fourier series. 

The complex Fourier series, defined over the interval (0,2p) for any 
arbitrary function, f(t), is given by: 

00 

f(t) ^ C^ exp(ni7Tt/p) (2.1) 

n=-«> 

where the complex coefficients, c^, are determined from: 

2p 

Cn = J f(t) exp(-ni7Tt/p)dt (2.2) 


The addition theorem may be stated as follows: 

If C„T and C ^ are the nth coefficients, respectively, of f (t) and 
n I n2 1 

f (t), then the nth coefficient, C„, off=f +f isC,+ C„«* 

2 n 1 2 n I n^i 


Proof: 


C 

n 


1 _ 

2p 


2p 


I 


f(t) exp(-niiTt/p)dt 



1 2p 

J exp(-niirt/p)dt 

0 

1 

C„ = — f f (t) exp(-niirt/p)dt 

2p J ^ 

0 

2p 

+ ]r~ r f (t) exp(-niiTt/p)dt 
dp j z 

0 

C - + C (2.3) 

n nl n2 

The shift theorem will be derived next. 

If the function is shifted by an amount a, the coefficients will be 
changed as follows: Let these coefficients be distinguished from the previous 

coefficients by primes. 

✓ 2p 

^n ^ Ip r exp(-niTTt/p)dt (2.4) 

^ -^0 

/ 1 

C = ~ f f(t-a) exp(-niira/p) exp(niiTa/p) exp(-niiTt/p)dt 

n 2p J 

0 

2p 

= — exp(-ni7ra/p) j f(t-a) exp(-ni (t-a)/p)dt 

Change variables: Let x = t-a dx = dt 

2p+a 

exp(-ni7ra/p) f f(x) exp(--niiTx/p)dx (2.5) 

n 2 p 

a 
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Consider the integral in Equation (2.5). 

2p+a 

_f f(x) exp(-niTTx/p)dx 
a 

2p+a 

+ f f(x) exp(-ni7TX/p)dx (2.6) 

2p 

The complex exponential repeats itself after every interval of 2p, and its 
value from 2p to 2p+a is the same as its value from 0 to a. If f(x) is a 
periodic function of period 2p, then its value from 2p to 2p+a will also be 
the same as from 0 to a (the shift theorem will only be applied to functions 
with this property). We may then write Equation (2.6) as: 


2p 

= j f(x) exp(-ni7Tx/p)dx 
■^a 


2p 


J* f(x) exp(-niTTx/p)dx + J' f(x) exp(-ni7rx/p)dx 


2p 

/ 

0 


f(x) exp(-ni7Tx/p)dx 


Substituting this result in the integral in Equation (2.5): 

1 2p 

C = TT— exp(-niiTa/p) C f(x) exp(-ni7rx/p)dx 
n ^p J 

0 


C = exp(-ni7ra/p) 
n 


_1 2p 

L f(x) exp(-ni7rx/p)d>J 


(2.7) 


C = exp(-ni7ra/p)C 
n n 


( 2 . 8 ) 
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Thus the nth coefficient for the shifted function can be obtained simply 
by multiplying the coefficient of the unshifted function by the phase factor 
exp(-ni7Ta/p). 

With this result the convolution theorem for complex Fourier series may 
be derived. 

The data given are assumed to be a convolution of resolved peaks with 
some "machine function" that causes a spreading and "smearing” of the peaks 
together. The complex coefficients, then, are the transform of this con- 
volution : 


C 


h 


n 



2p 




>1 

f(x)g(t-x)dx 

j 


exp(-ni7Tt/p)dt 


( 2 . 9 ) 


where g(t) is the machine function, and f(t) is the assumed form of the data 

before convolution. Let denote the coefficients of the data, the 

" f 

coefficients (transform) of the machine function, and the coefficients 
(transform) of the deconvolved data. 

Note that the continuous integral was used rather than a discrete sum 
as of the discrete series. This is because the actual convolution was 
assumed to be a continuous process. Actually, as long as we have sufficient 
data points, the results of both methods should be essentially the same. 

The integral permits convenient mathematical manipulation. 

Assuming sufficient continuity so that the order of integration may 
be interchanged: 


h 

C 

n 


2p 

/ 


r 




~ r f(x) 

2p 


2p 

1/ 

V 0 


g(t-x) exp(-ni'irt/p)dt 


dx 


J 
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Applying the shift theorem: 


h 

C = 
n 


c" = 

n 


2p I 1 

f f(x) exp(-niiTx/p) ~ f g(t) exp(-niTrt/p)dt dx 
Jq 2p J 0 I 

L y 


2p 


J f(x) exp(-niT:x/p)C^ dx 


= 

n 


2p 


C^2p ^ J' f(x) exp(-niTTx/p)dx 


dl = 


2pC^C^ 
n n 


( 2 . 10 ) 


The product of the coefficients of the original function and the machine 
function multiplied by the constant factor, 2p, gives the coefficients of the 
convolved function. When all that is given is the data and some machine 
function, and the separated peaks are desired. Equation (2.10) is usually 
written in the form: 


= 1 - — 

n 2p (>g 


( 2 . 11 ) 


The coefficients of the separated peaks may be obtained by dividing the 
coefficients of the data by the coefficients of the machine function and 
multiplying by the constant factor l/2p. Transforming to the time domain will 
give the resolved peaks. This process is known as "deconvolution". It should 
be apparent that there is no unique set of coefficients that will satisfy this 
equation. The coefficients of the machine functions used in this research 
die out rather quickly at the higher frequencies and are essentially zero over 
most of the spectrum. Assuming cjj is zero in these regions, we have the 
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indeterminate case and almost any spectral function given for these regions 

where is zero added to the unique determination c[]/c9 where is not zero 

will satisfy Equation (2.11). One important particular solution gives = 0 

when = 0. This is known as the "principal solution" as given by Bracewell 
1 ? 

and Roberts. 

Note here that a single isolated peak for the case of gas chromatography 
cannot be taken for the machine function as it is not the response of an im- 
pulse; it is the response of one of the original resolved peaks "smeared" by 
the machine function to give the data. It is not possible nor desirable to 
assume the original input were "spikes." The sidelobes quickly become formid- 
able as the peaks become sharper. (See Appendix A for a discussion of sidelobes.) 
In accordance with the hypothesis that the area under the peak is proportional 
to the amount of material present, the deconvolution needs to be carried only 
to the point of barely resolving the peaks, hence the assumption that the orig- 
inal distribution of material before "smearing" was in the form of barely 
resolved peaks. 

Of course, one realizes that the form of the convolution integral is 
obtained by considering the original input peaks to be broken up into infin- 
itesimal "impulses," each one of these giving rise to infinitesimal "machine 

functions" positioned in a continuous manner on the t-axis and that summing 

13 

all these together gives the data. 

A first guess at a good machine function, then, might be one of the same 
shape as an isolated peak, but much narrower; more specifically, if an isolated 
peak can be expressed in some functional form, g(t), then a machine function 
of the form g(at) would be tried, where a must be greater than one. This was 
tried with only moderate success. The best results were obtained by taking 
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this function and giving it more "skewness." The assertion that this function 
is a good machine function has no firm foundation in theory as will be discus- 
sed subsequently. First, the linearity and shift-invariance of the apparatus 
will be discussed. 

The flame ionization detector is linear over a wide range. Stated 
mathematically, one would say that for certain amounts of different gases run 
through the chromatograph singly and forming the distributions: 
f (t), f (t), f (t), etc., then differing amounts of these same substances 

1 2 3 

run through the machine together would form the composite distribution: 
af (t) + bf (t) + cf (t) + . . . 

12 3 

regardless of whether or not they were overlapping. The constants a,b,c, etc., 

are determined by the proportions of the substances in the sample. 

Even a cursory examination, however, reveals that the distributions are 

not shift-invariant. Being shift-invariant means that for an amount of gas 

coming out of the tube earlier or later than a given gas should have the same 

shape as that gas. To state this more explicitly, if the distributions are 

shift-invariant, then for an arbitrary distribution of functional form f(t), 

any other distribution could be expressed in the form c,f(t-C9), where c 

^ 1 

and c are constants. An examination of the data shows that the peaks seem 
2 

to be all of the same form but that the ones coming out earlier are narrower 
than those coming out later. Measurements show that the widths at half-max- 
imum become larger in a continuous manner as one goes from left to right in the 
data. This lack of shift-invariance will affect the accuracy of the results 
only slightly. This will be discussed in more detail subsequently. 

Taking the area under each peak as a measure of the amount of each 
material present, how this area would be changed, if any, under deconvolution 
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needs to be investigated. To look at this problem in its simplest form, free 
of all nonessentials, consider two isolated peaks of any arbitrary shape and 
size. The areas of each would be given by, respectively: 


2p 2p 



The ratio of their areas would be: 



Convolving each of these functions with some arbitrary machine function, g(t), 
would yield: 

2p 

hi(t) = r fi(t)g(x-t)dt 


for fi(t) and an expression of similar form for fz(t). 
areas after convolution would be: 



2p 2p 

t f fj(t)g(x-t)dt 


0 


'0 


dx 




2p r2p 

fz(t)g(x-t)dt I dx 






The ratio of the 


( 2 . 12 ) 


As they are both of similar form, only one expression will be considered. 
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2p f2p 




fi(t)g(x-t)dt 


dx = 


2p 


r 


2p 




r fi(t) r 

^0 Uo 


g(x-t)dx 


dt 




The expression on the right was obtained by interchanging the order of 
integration. The integral over g(x-t) will be a constant as it gives the area 
under g which will remain the same regardless of how much it is shifted along 
the axis by t. Let this integral be denoted by K. A similar result will be 
obtained for h^Ct) and Equation (2.12) may be written: 


2p 2p 2p 


r h,(t)dt 
■’o 

r f,(t)Kdt 
■^0 

fl(t)dt 

-'0 

2p 

2p 

2p 

r h2(t)dt 

r f 2 (t)Kdt 

■•o 

f2(t)dt 

-^0 

2p 

2p 


r hj(t)dt 

r fi(t)dt 


■*0 

''o 



2p 2p 


£ 


h2(t)dt 


£ 




(2.13) 


and the same ratio of areas as before convolution is obtained. Thus the 
same ratio of peak areas will be obtained under deconvolution. The case 
where the peaks are merged needs to be considered. It should be apparent 
after a little thought that if two or more peaks of unknown shape are 
merged there is no possible way to determine how much of the smeared together 
area to assign to each distribution. This situation is actually no different 
from the one where the peaks are well separated, but only more vividly illus- 
trates the problem. The deconvolution process requires closer examination. 
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From the convolution theorem, the coefficients of the deconvolved 


function are given by: 

= 1 - — ( 2 . 11 ) 

n 2p c9 


Suppose the areas are arbitrarily assigned under two smeared together peaks 
and the complete function written as the sum of two others: 
h(t) = hi(t) + h2(t) 

By the addition theorem the transform would be: 
n n, n^ 


h h 

where C and C are the transforms of hi(t) and h,(t), respectively, 

■^1 Hz 

Dividing this by the transform of the machine function, we get the coefficient 
of the deconvolved function. 


C 


f 

n 



1 _ 

2p 



1 _ 

2p 



+ 


1 _ 

2p 



C3 


n 


Taking the inverse transform and making use of the addition theorem, we get 
for the deconvolved function: 
f(t) = fi(t) + f2(t) 

where fi(t) and f 2 (t) are the deconvolutions of hi(t) and h 2 (t), respectively. 
Suppose the areas under the peaks were assigned such that two other quite 
different functions are obtained. 
h(t) = hgCt) + h4(t) 
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After deconvolution the result obtained is: 

f(t) = fait) + U(t) 

The deconvolved result is the sum of two quite different distributions. This 
is not a contradiction; only one result is obtained. For instance, if f^ 
and f^ were both zero on some region of the t-axis and f^ was negative a 
certain amount and f^ was positive by the same amount, then fg + f^ would 
equal zero in that region as would f^ + f2. giving the same result in each case. 
Of course, the areas under each of these distributions would in general be 
different, so that if one didn't know the shapes of the two original peaks 
before smearing, then there would be no way of determining the correct areas. 
However, if the correct shapes of the peaks were known, even if they were 
smeared together, then one could deconvolve them and know how much area to 
assign to each peak by examining the deconvolution of each peak separately. 

For this particular problem it is known within narrow limits what the 
peak shapes should be. All are reasonably uniform except for small variations 
in narrowness. If, then, a machine function can be found that gives an 
isolated peak a deconvolution with small negative regions, and if the decon- 
volved peaks all have similar form, then one may accept the results with a 
high degree of confidence. The requirement of small negative regions for the 
deconvolved peaks is if the deconvolution gives well separated peaks with no 
negative regions, then there will be no problem of cancellation of the area 
of adjacent peaks by these negative regions and hence there will be no problem 
in assigning the areas under the peaks. 

For the sake of completeness, a deconvolution of a linear, shift-invariant 
set of data will be included. 

Any number of distributions making up the data could be expressed in the 



following functional form: 


(2.14) 

Apply the convolution 


H(t) = ah(t-ti) + bh(t-t 2 ) + chCt-ta) + . . . 
where a, b, c, etc., ti, ta, etc., are constants, 
theorem, and employ the addition and shift theorems: 

p a exp(-niirti/p) exp(-niTrt 2 /p)cj^ + . . . 

" 2p ^ ^ 

" 2pc9 


c'^ = 

n 


/' > 
a exp(-niTTtyp) + b exp(-ni7rt2/p) + . . . 


1 c' 

1 n 
2p c9 


(2.15) 


where C" is the transform of h(t). Transforming back to the time domain would 


n 

yield the expression: 

F(g) = af(t-ti) + bf(t-t 2 ) + cf(t-t 3 ) + . . . (2.16) 

where f(t) is the deconvolution of h(t). 

So, if a machine function can be found that gives a minimum of negative 
regions for the deconvolution of an isolated peak, then the correct deconvo- 
lution of the data will give peaks all of the same shape, differing only in 
their heights. Note also that the location of the peaks will remain unchanged 
under deconvolution. 


19 



CHAPTER III 


DECONVOLUTION OF GAS CHROMATOGRAPHIC DATA 


The general method of performing deconvolution is as follows: One must 
initially obtain a suitable apparatus function. A transform is then made of 
this machine function and of the data. The transform of the data is divided 
by the transform of the machine function. Equation (2.11) is used for this. 


C* = — 
n 2p 

n 


( 2 . 11 ) 


These coefficients are then transformed back to the time domain to give the 
resolved peaks. The small imaginary part of these figures is due to calcula- 
tions roundoff and can be neglected. 

If the entire transform after division is transformed back to the time 
domain, the result, quite apparently, is nonsense. An inspection of the trans- 
forms reveals why this is the case (see Figure 3). Note from the graphs of 
the machine function and the data that the coefficients of the lower frequen- 
cies are the only appreciable ones (see Figure 4). They begin rather large 
at the lowest frequencies and fall off in a regular manner. It is apparent 
in nearly all cases that if this trailing off of the function were extrap- 
olated into the "noise" it would approach zero rather quickly, certainly in 
all cases before at most a hundred coefficients were taken. The coefficients 
falling to zero this quickly reveal that no information contained in our data 
was lost, getting at least two points per wavelength on the highest frequency 
present (see Appendix B where the sampling theorem is discussed further). By 
its relation to the continuous Fourier transform, it is apparent that the 
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I 

Fig, 3, Graph of the real parts of the complex numbers that result from divi- 
sion of the coefficients of the data by coefficients of the machine function. These 
coefficients are truncated at some point (see arrows) and transformed back to the time 
domain to give the deconvolution. 





Fig. 4. Graphs of coefficients of the data and two machine functions: (a) is 

graph of real part of complex coefficients of the machine function, (b) is real part of 
coefficients of data, and (c) is real part of coefficients of a gaussian function (not 


centered at zero) . 



envelope of the transforms of gaussians are gausslan also and would fall to 
zero rather quickly However, one notes by examining the transforms that 
a small residual "noise" extends the whole length of the spectrum. This is 
due to machine roundoff in the computer, noise and error in the original data, 
and digitizing of data. This noise seems to be completely random and as much 
positive as negative, so it is not necessary to raise or lower the spectral 
function to some "baseline." 

Digressing for a moment, the data and the manner in which it was recorded 
should be considered. Due to the mechanical and electrical properties, all 
recording instruments have a "smoothing" effect on the recorded data. Speaking 
in terms of the Fourier transform, one would say that the higher frequencies 
are filtered out. As discussed in the introduction, the area under the peaks 
is an accurate but not exact measure of the amount of material present. If 
the data were not smoothed and the recorder accurately measured the amount 
of material passing through an arbitrary cross section of the tube at every 
instant of time, the data would show statistical fluctuation due to the 
thermodynamic nature of the variables concerned. The smoothed data recorded 
is some sort of "average" of these statistical fluctuations and the error 
involved would in most cases be small compared to other inherent errors. 

Thus, for the purpose of taking the area under the peaks, the smoothed, or 
"bandlimited" form for the data and the deconvolved peaks will be sufficiently 
accurate for the purpose. If this is the case, the principal solution for 
deconvolution will be the correct solution. Practically, the presence of noise 
complicates the situation. The true spectrum of the machine function and its 
associated deconvolution buried in the noise will be unknown. This will 
necessitate truncation of the spectrum at some point before the true value of 
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becomes zero, 
n 

The reason the noise assumes such significance is because of the division 
of one frequency spectrum by another. One may accept with confidence the 
results in the region where both coefficients are appreciable. However, in 
the region of the noise, this division, due to the random nature of the noise, 
can give rise to enormous numbers. Note from the graph that the largest numbers 
are in this region. These results are completely meaningless and can quickly 
overwhelm the correct values if very many coefficients are taken as one extends 
into this region (see Figure 14). The results obtained when the dividing func- 
tion is small represents one of the largest problems encountered in practical 
deconvolution. It is quite apparent that these frequencies must be "cut off" 

before getting very far into this region. This is referred to as "window 

15 

closing" in the literature. 

Various methods of dealing with this problem have been suggested. Some 

merely accept the results as they are, taking it as being accurate enough for 

their purposes. Others, desiring more accurate results, employ the use of one 
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of the many "data windows" available (see Figure 5). ’ (See Appendix A 

for a discussion of data windows.) 

A Hamming window was applied in the frequency domain with only moderate 
success (see Figure 6). Most of the problem was due to widening of the trans- 
formed peaks. Such a strong deconvolution had to be taken in order to offset 
this widening that some distortion was introduced. There was no improvement 

over deconvolution with window not applied. Other workers in the field have 

1 8 

used data windows, however, with some measure of success. 

A data window, however, is an artificial device and its determinations 
only fortuitously resemble the original function because of its reduction of 
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A = 1 for t = 0 to T 



A= . 08+ . 46 ( l-cos27Tt/T) for 
t=0 to T 


Fig. 5. Four common data windows, 
(a) Rectangle. (b) Extended cosine bell‘. 


A= . 5 ( l~cos2‘iT5t/T) for T=0 to T/10 

and t=9T/10 to T 

A=1 for t=T/10 to 9T/10 



A=l-6(2t/T-l) ^+6l2t/T-l| ^for t= 
T/4 to 3T/4 

A=2 (1- I 2t/T-l I ) ^for t=0 to T/4 
and t - 3T/4 to T 

Data is taken over the interval (0,T) . 

(c) Hamming. (d) Parzen. 



Fig. 6. Deconvolution with Hamming window applied compared to straightforward 
deconvolution: (a) is data of Run 5, (b) is deconvolution with Hamming window 

applied, (c) is deconvolution with function continuation, and (d) is straightforward 
deconvolution . 



the sidelobes. The peak shape is also usually changed slightly. To achieve 
any substantial improvement in deconvolution, more precise information about 
each particular problem is needed. Specifically, if one knew exactly what the 
coefficients were supposed to be beyond the point of cutoff, the perfect 
deconvolution should be obtained with no sidelobes. A measurement of the wave- 
length of the sidelobes showed that it corresponds closely with the wavelength 
of the cutoff frequency. This suggested that the sidelobes and other distortion 
were possibly caused by "cutting off" the frequencies at this point. However, 
this turned out to be only partially correct. Most of the error was found to 
stem from a poor choice of the machine function. This will be discussed later. 
For the following discussion, the assumption will be that the deconvolution 
was performed with some "ideal" machine function. This was realized in 
practice, to a large extent, by later deconvolutions. 

Even with a good machine function, truncation of the frequencies will 
cause error; there will be sidelobes. The main thrust of the ensuing argument 
is that if a function could be found which resembled the deconvolved function 
closely enough, such that the coefficients matched up with the coefficients 
of the original deconvolved function, especially in the region close to the 
point of cutoff, then the coefficients should be close in numerical value to 
what the coefficients of the original function would have been had they not 
been cut off. To test this hypothesis, an artificial function was built from 
gaussian peaks of the same height, location, and width at half-maximum as the 
deconvolved function (see Figures 7 and 8). The transform was taken and the 
frequencies were cut off at the same point as the original function to see what 
effect this would have. The results were very encouraging. The sidelobes 
appeared at the same positions as the original deconvolution and had almost 
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the same amplitude. This indicated that most of the error remaining in 
the deconvolved data was due to the act of "cutting off" the frequencies. 
Coefficients of the artificial function were then substituted into the 
original transform from the point of cutoff on. The transform was then a 
hybrid with all of its original coefficients up to the point of cutoff, and 
with coefficients from that point up to a hundred coefficients "grafted" from 
the artificial function. The rest were set to zero as they originated in 
noise anyway. Transforming back to the time domain to observe the results: The 
results were very good for Run 4 (see Figures 7 and 9). The amplitude of the 
si delobes were reduced by more than half. The results from Run 5 were not 
quite as impressive, only slight improvement being gained. However, the 
coefficients after forty-six terms do not contribute much "ripple" anyway 
(see Figures 8 and 10). A little more work on the machine function is probably 
needed here. 

More work needs to be done in obtaining a good artificial function also. 
The matchup between the coefficients of the artificial function and the actual 
function was crude (see Figure 11). This suggests a more fundamental dif- 
ference between the two than a comparison of the data reveals. Encouragingly, 
though, the match was as good in the region close to the cutoff point as at 
any other region. The actual peaks were different from guassians in that they 
were narrower close to the base. A promising field of research would be the 
investigation of the efficacy of various peak shapes for artificial functions 
in deconvolution. One item that needs to be pointed out here is that the 
artificial function should resemble the original deconvolution after cutting 
off its frequency spectrum, not before. This would involve some adjusting 
as the solution is converged to. At any rate, the match worked well enough 
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Fig. 9. An enlarged view of the first three peaks of 
Run 4: . (a) is data of Run 4, (b) is artificial function, 

(c) is deconvolution with function continuation / and (d) is 
straightforv/ard deconvolution. 
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Fig. 10. An enlarged view of the first four peaks of 
Run 5: (a) is data of Run 5, (b) is artificiaJ function, 

(c) is deconvolution with function continuation, and (d) is 
straightforward deconvolution. 
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Fig. 11. Coefficients of the artificial and deconvolved function. 
"HOUT" contains the coefficients of the deconvolved function and "AROUT" 
contains the coefficients of the artificial function. 
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in practice to considerably improve Run 4. 

Also, there is no reason why substitution of coefficients should begin 
at the point of cutoff. More and more coefficients from the artificial function 
may be substituted until a more desirable result is obtained (see Figures 12 
and 13). However, one must keep in mind in this case that the final peaks will 
take on more of the character of the artificial function. 

Finding a good machine function will probably be the largest problem 
encountered. A gaussian for the machine function was easiest to implement, but 
left curious "shoulders" on the right side of some of the peaks, giving the 
impression that there were additional peaks that could be revealed by taking 
a stronger deconvolution. This was highly unlikely, however, as there was no 
other evidence for this. It was suspected that the shoulders were due to not 
considering the skewness of the peaks. 

In the next attempt, several isolated peaks in Runs 1 and 2 were tried 
as machine functions. Narrowing the peaks was easy to implement in a computer 
program by selecting every other, or every third data point, or some other 
interval. Interpolation between the points was introduced later to obtain 
a finer gradation. The amount of resolution desired in the final deconvolution 
is obtained by adjusting the narrowness of the machine function. 

However, only small improvement over a guassian was gained by using 

isolated peaks (see Figure 14). It was apparent that a function with more 

"skewness" was needed. There are several functional forms, all referred to 

generically as "skewed gaussians" that have been used successfully by other 

1 9 

workers in the field. The method used with success in this research was to 
choose an isolated peak, which was close to the desired shape anyway, and give 
it an additional "skewness" (see Figures 7 and 8). This was accomplished 
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Fig. 12. Comparison of two different cases of func- 
tion continuation for Run 4: (a) is data of Run A, (b) is 

deconvolution with function continuation (20 to 100 coeffi- 
cients substituted from artificial function) , (c) is decon- 

volution with function continuation (35 to 100 coefficients 
substituted), and (d) is straightforward deconvolution. 


37 



Fig. 13. Comparison of two different cases of func- 
tion continuation for Run 5: (a) is data of Run 5, (b) is 

deconvolution with function continuation (37 to 100 coeffi- 
cients substituted from the artificial function) , (c) is 

deconvolution with function continuation (47 to 100 coeffi- 
cients substituted) , and (d) is straightforward deconvolu- 
tion . 
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Fig. 14. Deconvolution with *a gaussian and an isolated peak as machine 
function. Also distortion in deconvolution introduced when too many coefficients 
are taken: (a) is data of Run 4 , (b) is deconvolution with coefficients truncated at 

60 coefficients (distortion is considerable, (c) is deconvolution with 5th peak in 
Run 1 as machine function, and (d) is deconvolution with a gaussian as machine function 




as follows: Taking an origin at the center of the peak, the function was 

compressed behind the origin and expanded in front of it, the compression 
and expansion being accomplished by only one factor (see Figure 16). A com- 
puter program was written to implement this and also to adjust the narrowness 
of the peak. By using a convenient feature of Xerox extended FORTRAN, the 
INPUT statement, three parameters were input after the program was run which 
seemed to be the only ones necessary for good deconvolution. The first input 
controls the narrowness, the second adjusts the skewness, and the third shifts 
the function either left or right on the t-axis. Various combinations of 
these parameters are tried until a good deconvolution is obtained. What 
constitutes a "good" deconvolution is more fully discussed in Chapter II. If 
the si delobes are comparable in size to that obtained from the artificial 
function when its coefficients are truncated at the same point, then one has 
good reason to believe that a limit has been reached beyond which one cannot 
go without using function continuation. 

Also, one should keep in mind that a good deconvolution may not be unique. 
To reiterate, a "good" deconvolution for the purpose of gas chromatography is 
one in which the area associated with the deconvolution of an isolated peak 
is coalesced into one well defined "peak" and the negative regions are small 
so that there will be no problem in assigning the areas under the peaks. 

There may be a number of machine functions and their associated deconvolutions 
that satisfy these criteria. 

An important consideration is that a good deconvolution seems to depend 
critically on small irregularities in the machine function (see Figures 15 
and 17). Note from Figure 15 that when the machine function that drops to 
zero in the skewed portion is used, the graph of the deconvolution has larger 
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Fig. 15. Comparison of a smooth machine function to 
one with a sharp discontinuity. On left is machine function 
as program MACHFUN outputs it. Note the dropoff in the tail 
On right is corrected machine function. 



Fig. 16. Comparison of a skewed with an unskewed 
machine function. On right is 5th peak in Run 1. On left 
is same peak with "skewing." This was machine function 
used for Run 4. c = .4 for this peak. 



Fig. 17. Effects of irregular machine function on 
deconvolution: (a) is data of Run 5, (b) is deconvolution 

with irregular machine function (left peak in Fig. 15)/ and 
(c) is deconvolution with corrected machine function (right 
peak in Fig. 15) . 



si delobes and shows more distortion than the case where the machine function 
is continued in an approximately linear fashion to zero. 

To obtain a good deconvolution involves a certain amount of trial and 
error. An inspection of the outputs at each stage is usually necessary. 

Once good deconvolved peaks were obtained, the areas under the peaks 
were acquired by numerical integration. Because of the sidelobes there were 
usually small negative regions on both sides of the peaks. The areas under 
the peaks were obtained by integrating between these negative regions as no 
other methods were known for accurately treating this particular situation. 
The suspicion is that if the sidelobes did not exist the peaks would probably 
be si ightly-broader at the base and that the true area would be slightly 
greater than the values determined. 
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CHAPTER IV 


ERROR ANALYSIS AND CONCLUSION 


Error Analysl s 


The trapezoidal rule was used to determine the areas under the peaks. 

^1 1 1 

2 ^ 2 3 2 ^-1 


Area = 


At 


(A.l) 


The data points were so close together that any difference between this and 
Simpson's rule would be so small that it could safely be neglected. Also, 
the endpoints were either very small or zero, so that an excellent approxima- 
tion can be obtained by simply summing the values of the data points. The true 

areas would be obtained by multiplication by the proper proportionality 
20 

constant. 


Area = prop, const. 


r 

Xi + Xi + X3 + 
L 



(A. 2) 


The results of any particular run are usually normalized; that is, the 

percentage composition is calculated by measuring the area of each peak and 
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dividing the individual areas by the total area, e.g.. 


% B = of B . ioq 

Total Area 


(A.3) 


The normalized calculations for all four runs are given in Table 2. The 
agreement of the percentage ratio of each peak among runs is relatively good. 


even for the deconvolved results. 


Areas of compounds are not directly proportional to the percentage compo- 
sition, i.e., different compounds have different detector responses; therefore, 
it is necessary to determine correction factors. Once determined, these cor- 
rection factors can be used to calculate the percentage composition.^^ 

In an attempt to determine some idea of the inherent errors in the taking 
of gas chromatographic data, ratios were taken of the peaks in Run 1 and com- 
pared to the same ratios in , Run 2 to see if there were significant differences. 
Only Runs 1 and 2 were compared because the peaks were well separated here. 
Although most compared quite well, one difference of 5 percent was observed. 
This means that it is difficult to get a very exact error analysis. 

Other methods of determining the amounts of material present from over- 
lapped peaks are triangulation, peak height measurement, and the dropping of 
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a perpendicular at the lowest point of the valley between two peaks. 

Obviously, none of these methods could apply to the first two peaks of Run 5 

as the peaks are almost completely merged. This leaves deconvolution as the 
only method available for treating problems of this type. In Run 4 the first 

two peaks are merged to such an extent that none of the above methods are 

recommended.^'^ The perpendi cular drop method will be investigated to get 
some idea of the error involved. From Table 2 the ratio of Peaks 1 and 2 in 
Run 4 results in a value of 2.47. Deconvolution gives a value of 1.57. The 
ratio of the same peaks in Run 1 where the peaks are well separated is 1.64. 

The ratio of the deconvolved peaks differs from the ratio of the same peaks 
in Run 1 by 4 percent. The ratio from the perpendicular drop method differs 
from the ratio in Run 1 by 51 percent. This is much larger than one would 
expect from experimental error. In such cases one would prefer the 
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TABLE 2 


AREAS UNDER THE PEAKS 
(Normalized) 





Percentage Composition 


Sample Mixture 

Run 1 

Run 2 

Run 4 Run 5 

After Decon. 

Run 4 Run 5 

Before Decon. 


Compound A 

13.010 

12.768 

13.372 

13.418 

9.814 

{33.840 

Compound B 

21.371 

21.428 

20.960 

20.827 

24.219 


Compound C 

33.856 

33.470 

33.082 

33.115 

32.125 

31.950 

Compound D 

14.223 

14.656 

14.714 

14.531 

15.380 

15.346 

Compound E 

17.541 

17.678 

17.871 

18.108 

18.462 

18.865 


deconvolved value if it were available. 


Discussion of Results and Conclusion 

One of the goals of this project was to develop methods for obtaining 
a good deconvolution with a minimum of time and effort. This was achieved 
in large measure by developing computer programs to handle all aspects of 
the problem and using INPUT statements to most conveniently vary the impor- 
tant parameters. 

Striving for generality was another goal. Working in the timesharing 
mode and using FORTRAN exclusively will not be the methods everyone will use. 
However, the programs will not be too difficult to adapt to the minicomputer 
or other data handling machine the researcher may have. The methods used 
here will not have general applicability. They would be most appropriately 
applied to data where the deconvolved result would consist of a number of 
quite distinct "peaks." However, in cases where the machine function is known 
exactly, the concept of function continuation could probably be applied in 
some sort of "iterative" process that converges to the correct result. Perhaps 
the sidelobes could be identified to give a starting point. 

The concept of function continuation in the frequency domain provides 
a twofold advantage. Not only does it serve to improve the deconvolution 
considerably, but, by examining the truncated artificial function, a basis for 
determining how close one is to the attainable limit in a straightforward 
deconvolution is obtained. 

More work needs to be done in making good artificial functions as evidenced 
by the crude matchup of coefficients. That the peak shapes are not gaussian 
is obvious when the artificial function is overlaid over the deconvolved result. 



One needs to find some distributions that can be expressed algebraically 
that closely resemble the peak shape. Alternatively, the shape of the decon- 
volved peaks may be changed slightly by changing the shape of the machine 
function. 

The shape of the machine function used in this investigation was obtained 
by taking the top of a given peak as the origin, compressing the t-axis behind 
this origin and expanding the t-axis in front of it. The amount of compression 
and expansion was varied by inputting a value for the variable "c" in the 
computer program. This was expressed in the program as (x+c*x), which is 
algebraically the same as (c+l)x; c+1 is the linear factor of expansion 
and compression. It took the value 1.4 for the deconvolution of Run 4 and 
1.36 for Run 5. It is easy to see how this might be generalized to give the 
t-axis any variable amount of expansion and compression, and hence give the 
machine function and shape desired, by using a general polynomial in place 
of this linear factor. 

a+bx+cx^+dx^+ ... 

The constants a, b, c, etc. are varied to give the results desired. A 
different polynomial may be used for expansion from the one used for compres- 
sion. Note that the expression above is a special case of this general 
polynomial with b=c+l and a=c=d= ... =0. 

Some of the programs are inefficient as far as computer time is concerned. 
No attempt has been made to make them more efficient nor to write them in 
more polished form. In particular, if it was all to be done over, the pro- 
gram for the machine function would be rewritten using a guassian rather than 
an isolated peak as a beginning function; thus eliminating the need for 
interpolating between the points. 


48 


APPENDIX A 


DATA WINDOWS 

Data windows will be discussed only qualitatively here. For a more 

rigorous mathematical discussion, one is referred to the references 

The frequency and time transforms are very similar in form, differing 

only by sign on the exponential term. The results obtained are nearly always 
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of the same form regardless of in which domain we begin initially. Thus, 

a data window will have essentially the same effect regardless of the domain 

in which it is applied. It is better interpreted when applied in the time 

domain so attention will initially be restricted to this. 

A knowledge of Fourier transforms shows that the sharper points of a 

function require high frequencies to represent it, an abrupt discontinuity 

taking frequencies ranging to infinity. The Fourier series is periodic so 

that if the function and its slope are not the same at its endpoints there 

will be an abrupt discontinuity and the frequencies will range to the highest 

ones present as the repeating series tried to join itself over the endpoint 

values. These frequencies usually group around the major frequency components 

27 

making up the function and is referred to as "leakage" in the literature. 

(See Figure 18.) 

To more faithfully represent the frequency spectrum over the real, infinite 
time domain would require a data window over a longer interval of time; the 
wider the data window, the better the results. However, this is not often 
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(b) Function and its slope the same at the endpoints. Transform is two 


symmetrical spikes. 


Fig. 18. Cosine function and associated transform with data taken over two 


different intervals to illustrate leakage. 




possible to realize in practice. Also if, by chance, the function and its 
first derivative happened to be the same at both endpoints, these additional 
frequency components would not be introduced. This happy situation seldom 
occurs naturally, however. One may, in a sense, bring this about, however, 
by extending the function at its endpoints in some smooth curve such that the 
function and its first derivatives are smooth and continuous across this 
region. Another way of reducing these additional frequency components is 
by using a "data window." 

Actually, the operation of taking data over a finite range is referred to 
as applying a "rectangular" data window as it can be viewed as the multiplica- 
tion of a rectangular function of unit height times the actual time function 
over the real, infinite time domain. 

Data windows are functions that decrease to zero outside the range over 
which data is taken and are multiplied by the real, infinite time function. 
They all usually have a value of unit in the center and fall off in a regular 
manner to small values or zero at the endpoints. In all cases their purpose 
is to make the function and its first derivatives smooth and continuous at 
the endpoints. 

Note that the frequency spectrum introduced by a sharp discontinuity 
usually expresses itself in the form of "sidelobes" on both sides of the 
function. (See Figure 19.) 

Fortunately, the data in this research is not beset with the problem 
of abrupt discontinuities. The data and its slope go smoothly to zero at 
both ends of the data. 

In the frequency domain, however, the researcher is confronted with 
this abrupt behavior as the frequency spectrum is cut off at some point. 


(a) Rectangle function, II (x). 



(b) Transform of rectangle function, 
sine X. Note the sidelobes. 


Fig. 19. Rectangle function and its transform, 

sine X. 
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As might be expected, si delobes are obtained with the transformation to the 
time domain. A data window has the effect of bringing the frequency function 
and its slope relatively smoothly to zero or small values, and, hence, to reduce 
the si delobes. As the data window tends to narrow the function somewhat, the 
effect in the transform domain is to broaden the peaks slightly. This broad- 
ening, which varies with the type of data window used, is one of the worst dis- 
advantages of these windows. Also the shape of the peaks may be distorted 
si ightly. 

The above statements, apparent for a single distribution, may not be so 
obvious when the data consist of a number of peaks at differing locations. 

But if one considers that the transform of a sum is the same as the sum of the 
transforms (addition theorem) and that a shift in location means only multi- 
plication of the transform by a "phase factor" of unit amplitude, then it is 
easier to see how all the peaks will be affected in the same manner. 
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APPENDIX B 


THE DISCRETE FOURIER TRANSFORM 

The complex Fourier series representation of a function defined over the 
interval (0,2p) is given by: 


f(t) = I C„ exp(niirt/p) 


( 2 . 1 ) 


n = -00 


where the complex coefficients, C^, are given by: 


2p 


Cn = ^ f f{t) exp(-n1irt/p)dt 


( 2 . 2 ) 


If only discrete data points are given the above integral can be 
numerically approximated by a summation. Suppose there are N equally spaced 
data points, Xj^: 


t = kAt 
1 


C = 
n 


N-1 

y X exp 
NAt k=0 k 


NAt = 2p 
'^niirkAt ^ 


n4^ 


2 


At 


At - dt 


C = — 
n N 


N-1 

I X, exp(-ni2TTk/N) 
k=0 


(B.l) 


If the data are band-limited, that is, the frequency spectrum does not 
extend past a certain point, and if the sample point interval is smaller than 
a half-wavelength of the highest frequency present, then no information will 
be lost in taking the summation rather than the continuous integral in 
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determining the coefficients. If the above is true, then more than N 
terms need not be included in order to completely restore the original 
function f(t). 


f(t) = 


2 

I 

n = 


C exp(ni7Tt/p) 
n 
N 

2 


(B.2) 


If discrete points are desired rather than a continuous function, this 
may be written as: 


X 

k 



'^iirkAt'^ 


exp 




N ^ 

2 J 


N - 1 

\ C exp(ni2TTk/N) 


(B.3) 


By considering the periodicity of the exponent and by noting that the deri- 
vation of from -N/2 to -1 is the same as from N/2 to N-1 , it is apparent 
that this expression may be written: 

N-1 

X. = I C exp(ni2TTk/N) (B.4) 

^ n=0 


The expression in this form is called the "inverse discrete Fourier transform" 
and its companion: 


C = 1 y X exp(-ni2TTk/N) 
n N k 


(B.5) 
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the "discrete Fourier transform." 

The equations expressed in this symmetric form make it possible to write 
only one computer program to handle both computations. One has only to change 
the sign of the exponent and to consider the factor 1/N when doing the inverse 
transform. 

The fast Fourier transform, usually written FFT, is only a computational 
algorithm to shorten the length of time over conventional methods in comput- 

P7 pq “30 

ing the discrete Fourier transform.^' 
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APPENDIX C 


COMPUTER PROGRAMS WRITTEN TO PERFORM 
FFT DECONVOLUTION 

Due to the enormous volume of calculation involved, computer programs 
were written to handle almost all aspects of the program. 

The computer used was the Xerox Sigma 9 which operates under Control 
Program V (CPV version EOO). Most computing was done on the cathode-ray 
terminals (time-sharing mode). The Tektronix 4010-1 graphics terminal and 

a 

a Calcomp 565 drum plotter were used to plot all graphs. All programs 
were written in FORTRAN IV with the exception of two Xerox extended FORTRAN 
options, the INPUT and OUTPUT statements. They serve as free form read and 
write, respectively. 

No subroutines were used with these programs. In the time-sharing mode 
it was more convenient to use the SET statements of the terminal executive 
language, especially in the earlier stages of the research when it was 
necessary to examine the output of each program. It should not be difficult 
for the programmer to adapt the programs to any system, however, with slight 
modi fi cations . 

FFT deconvolution was implemented with these programs as follows: 

First, a machine function was generated using either MACHFUN or TESTFUN. 
Program TESTFUN generated data points from a gaussian distribution. The 
height and width of the gaussian were varied by inputting two variables via 
the INPUT statement. MACHFUN first read in data points from the 5th peak in 
Run 1. Then by interpolating between the points the narrowness was varied in 

Tommy Dearmond, personal correspondence. 
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a continuous manner. The skewness was varied with one parameter by compres- 
sing the function behind the peak and expanding the function in front of it. 

The third parameter enabled the function to be shifted along the t-axis. This 
has the effect of shifting the deconvolved output. This was necessary in order 
that the program that builds the artificial function would not have "split" 
peaks to contend with. 

Next, a forward transform was taken of the machine function using program 
31 

FFPEAK to get its coefficients. Also, the forward transform was taken of the 
data using FFPEAK. 

These coefficients were read in program FGH and the coefficients of the 
data were divided, one by one, by the coefficients of the machine function. 

When applying a data window, this division is multiplied by a factor which 
depends on the data window applied. See program FGHAM to see how a Hamming 
window was applied to Run 5. 

The coefficients obtained by division are read into FFPEAK and an inverse 
transform is taken. The user is allowed to select the number of coefficients 
wanted to go into the transform via an INPUT statement. The user can thus 
cut off the frequency spectrum at any point and note which gives the best 
deconvolution. 

The real part of these numbers is read by program CALCOM and plotted. 

(If the imaginary part of these numbers is very large, one should suspect 
an error.) This is usually plotted along with the original data, and usually 
the hybrid function and artificial function, in order to better compare them. 

Program SUB is used to construct the hybrid function. Coefficients from 
the artificial function are substituted in beginning at any desired point. 

To use this program, input the number of coefficients of the original function 
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that is wanted and the rest will be substituted from the artificial function. 
An inverse transform of this is taken using FFPEAK. The result is usually 
plotted along with the straightforward deconvolution in order to better 
compare them. 

ARTFN is the program that builds the artificial function. It reads the 
deconvolved function and determines the location, height, and full width at 
half-maximum of each peak. An artificial function is then constructed with 
gaussian distributions having these parameters. This function may be plotted 
directly using CALCOM. The forward transform is then taken using FFPEAK 
to get the coefficients. 

In concluding, there is one point that needs clarifying. In some of the 
programs one, loosely, is asked to input the number of coefficients of that 
number of terms of the lowest frequencies. With the exception of the first 
coefficient, there are two coefficients that correspond to each wavelength. 
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NAME - TESTPUN 


C this generates data points CF a f^AUSSlAN 

C DlSTRI3Jri9N. there ARe TWO INPUTS# THE FIRST 

c oeteRumes the heIQht of the Peak and the second 

C DETERHI^^E3 THE WIDTH* 

C 

DIHEWSI3N Y<256> 

INPUT S,T 
DO 1 

Y< I )«3*3EXP<-(I-129. )»¥2»T) 

1 CONTINUE 
wRlTij 3*2 >Y 

2 F0R"1AT( 10F5*0I 
STOP 

end 
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name ■ ' 1 A*W^U‘J 


C QIveN ASy INITIAL SET OF Oaia POiMS OF i PROSPECTIVE 
C MACHINE pjnCtION, THIS PROORaH, BY USlNn ThREE INPUT 
C PARAMETERS/ GIVES THE FUNCTION a CERTAIN AMOUNT OF 

C 9KEHNE33, n-ARROws OR UROAOENS THE FUNCTION/ ANO SHIFTS 

C THf FUNCTI 9n later ally. In THAT ORDER. 

C 

dimension RUnIN( 102A )/MF« 256 )/»RUmo2A) 

DATA MF/0/ 

REA^ ( 35/ 1 » RjNIN 

1 FORMATj 10F5*0) 

INPUT C,SC,IBeQIN 

c 

C T^£ F<JLL3/*1N3 9It OF CODg; CCMPReSSeS THe Y.aXIS BeHINO 

C the P£A<. POINT 433 IS AT TH£ TOP OF THe PEAK. 

C 

DO 2 1-1/1024-433 
RU( Xf433»-RUNINl I ♦433* 

X-I 

K-INT( X-C»X» 

run IN I I ^433 )-| X-C*X.K )*IPU(K4*3A) 

A-RUK + 433) )*RU(K+4 33> 

2 CONTINUE 
C 

C this bit 3F code expands the X.AXIS IN FRONT OF THE PEAK, 

c 

DO 3 I-*»32/l/-l 
^-433-1 

K»433-1NT<X+C*X> 

IF ( <,uT,HRUNlN( I )-3 jGO TO 3 

RUNIN( I »-BUNIN(K»-{ ( X+C • X ) - I N T ( X*C *X ) \* 

#( HUNIN( <) -RUNIN(K-1 ► » 

3 continue 
C 

c the f0i-L34INQ bit Op cODE NaHROwB AND SH*pTS THE FUNCTION. 

C IT can obtain Data points for any Increment By interpolating 
c between the points. 
c 

UO 3 1-1/ INTt 1024/SC ) 

RI-I 

K-RWSC 

5 MF( IflBEQTN)-RUNIN(K) + (Rl*SC-K)*(RUNIS(K + H-RUNlN(K) J 
mRITE( 36/41MF 

4 FflRMATdOFS.O) 

STO» 

END 
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non o no n 


N*MC - PrPsA< 


c 

C THIS COHPUTES THE FBRWaKD OR INvfRSE FOURJER 

C tRaNSFORH POR AMY S£T CF OISCRETt OaTA POINTS WHICH 18 AN 

C integral »9wfr of two* input - 1 . for a forward transforh 
C and 1 , FOR AN inverse transform. IF . IS INPUTTED THE 

C PRoGRaI will read only real DATA POINTS oF A R£T OF ANY 

c length, IF 1 , IS inputted the program will read only 
C PRoPeRLY F9RHatTeO complex numbers* any arbitrary 

C SET OF DATA, BE SURE TO DIMENSION THE PROPfcR ARRAYS 

C TO the SHaLLfST POWgR OF TWO THAT WILL CONTAIN THE DATA. 

c tHe Rest op the values in the arRay will be set to zero 
C TO give an flwePALL data field of length an integral power 

C OF TWO, SET NSTAQE TO THIS INTEGRAL POW£R OF TWO, 

c aLso^ for thf inverse transform you will Need to input 
C T^e number of coefficients you WISH; tHe Re 8 T WILL BE S£T 
C to ZERO. 

E note: although the pResent ARRAvi are Dimensioned to hold 

C lOgA data points, the PROQRaM IS SET TO READ ONLY 356 * 


THIS first bit of cODE IS TO FI^U FHE ARRAYS. 

REAL*B 3 ATAC 102 A)*CeEF( 2 ,l 024 j 
REALMS FLTN,pHlaN,SlQN,FL IN2J*TEMP 
DOUSLE cOmPLEX W/X<g,1024) 

DATA DATa/0/ 

INTEGER R 
N8TA3E«lo 

SIGN--1 . 

N«2wi^N3r age 

IF ( 513N*ST«0* 1QO TO 33 

R£AD( 11* 77»END«6S> (OATa ( K ) *X«38 5, 640) 

77 FORMATnOFS.O) 

66 00 22 I-l.N 

22 X(li I ).DCMPLX(OATA( I >*C* > 

GO TO 55 
33 INPUT M 

H£AD( 13. 9S)C0EF 
99 FORMAT ( D2l •16«2X,D2 i .15 ) 

DO 9 1 -M+liN+l-M 
CflEF( 1 , n» 0 . 

9 COEFcg.D.O. 

DO 44 I* 1 *N 

44 X(l* 1 ) -DCmPLX(CDEF(i, I )*C0EF(2, I | ) 

S5 CONTINJE 

actual CALCUI ATION STARTS HERE. 

N2-N/2 

fltn«n 

PH 1 2S«6.8R3t85307179586/FLTN 
DO 3 J«liNSTAQE 
N2 JaN/ ( 2 a«J ) 

NR-N2J 
NI« ( 2** J ) /2 
00 2 I«1.NZ 
Ts? ,1. f T > *N? J 
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FLI'42JalN?J 

TEriP-r_,lN?J*FHl2N*8lQN 
W«DC1PLX<DC0S<TEMP)#0SIN(TEMPJ » 
oa 2 P>1*NR 

lSUBl.^tlN2J*2 

ISJ92«I3uni4NsJ 

ISU33*ISUB-t-N2 

X(2« I3J3 laXCl# I8UB1 )4W»X( 1> ISlf^BZ } 
X(2>I3JBai«X(t4 ISUBl )«U«X t 1# ISUB2) 

2 CONTiNJt 
DO 3 RalfN 

3 X(1«9)«X(?«R» 

1F( SIQ-J.GT.O. )Q6 TO 6 
DO *» 

♦ X(2«R)*X(i«R)/FLTN 
5 WRITE) 12«88) (X(2«J)iJ>l«N) 

81 FORI&ri 32i*l5«2X«D2l*l5) 

STOP 

END 
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name • 


C this PR33^AH divides THE Ct5EFP IC XE^TS OF THE TRANSFORM 
C OF the 3WEV DATA POINTS BY THE COEFFIClrNTS OF THE 

c machine fjvction Term by term to give the coefficients 

C OF the deconvolved FUNCTION, 

c 

c 

DOUBLE Complex F(102AI,G(10EAI,H(102ai 
REA3(8> 1 

1 FORMATiD2t.i5,Ey,OEi.i5) 

2 READOf3)H 

3 format ( D2i , ig,2X«D2l .I 5 ) 

4 DO 5 Iat«l024 

5 F(H.H(I»/G(I) 

MRI TEl 10»^ )F 

^ FORMAT ( D?4 .15, 2X,D2l .15 ) 

STOP 

END 
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n n 


name - C4L = 9«1 


C 

C THIS 19 the PR8QRAH THAT INITIATES TH£ PLOTTING ROUTINES, 

C IT WILL RcaO A9 HANY AS FOUR DIFFERENT OaTA FILES AND 
C PLOT THEH ALL ON ONE QRAPHt THERE ARE THREE INPUTS, FflR 
C the first, input EITHER A OR A »0», A »2* IF YOU 

C WANT the SjTPuT ON !«£ TEKTRONIX SCREEN, OR *01 IF YOu WANT 
C output on the CaLCOMP plotter* TiE SECOND INPUT IS A SCALING 
c Factor, the proqrah is ocsiqneo to graph all numerical 
C values in inches, all numerical values TIHeS the scaling 
C factor will give AS TH£ OUTPUT tHIS PRODUCT IN INCHES, THe 
C third input is the amount in DatA-POINT units that some of 
c the qRaRhs aR£ shifted, this is Necessary aeCause a shift 
c IN the machine Function cause® a shift in the output, this 
C must Be adjusted in order to CoMHaRc it with the ORIQINaL 
C data* 
c 

DIMEV3I9N X(lo24)#Y<lo?'*)*YY(i02MUZfio2*> 

RE*L*8 EOVERV 

Common /model/ itek 
INPUT ITEKaFCTR 
CALL FaCT»R(FCTRJ 
A EOVERV-0. 

R£A3(15,55)y 
55 format ( 321 «15 ) 

READ{ 17M11 IZ 
111 FORNATC 121 .15 J 

REA3( i6*6&#END>77 ) VY 
66 FORMAT! 10FS*0) 

77 DO 11 1-971*1026 
11 YYCn.O* 
iNPjT M 
00 B 1-1,1026 
X| I J-EOVERV+,5 
YY< I ) -YY! T l«,003-f5« 

Y( I ) -Y( I )¥*026+1* 

2( I ) -z( I )-*026*2,S 
100 1F< YY(1 >,nT,10* )YY(H-9» 

1F< YY( I > ,LT,0, »YY( I )-.00l 
IF(Y(I),GT,10,)Y(I)-IC. 
lF(r(l),LT»0,)V<I »-.00l 
EOVERV-1/60, 

B continue 

Call pl9T8(i,i»d 

CALL A<ISt,5*0*«lHT,l,26.*0,*6,,0*,l . ,6 hF6,1 ) 

CALu 4Xl6{ ,S>0,#1HG, 1,10,,90«*1,*0,#1 »*6 hF6,1 ) 

CALL PLar(X(l)*Y<M),3) 

DO 9B <-Z,1026-H 
68 CALL 9L0T(X(KI*Y(K-*M-1 )*2) 

00 22 <-1O25-m*1026 
22 CALL PL0T|X(K)*Y!<+M-1026)*2) 

CALL PL9T(X(i I* 2{i)i3) 

DO 99 <-2,1026-H 
99 CALL pL0T(X(K)»Z(K4M-1)#2) 

00 95 <-1020-M* 1026 
95 CALL PLOT(XCKI*2(K¥M-102a1 «2> 

CALL PLOTtXd }*YY(1 »f 3 I 
00 66 <-2,1026 
CALL 3l9T(X<K J,YY«X 1,;^ ) 
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C CALL LI'^£«X,Y,1024, 1, A, .2,0) 

C Call LINtrcX, yYa1024,1#4, .2/0> 

CALL STflPPLOT 
STOP 
END 
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noon 


name » ARirN 


c 

C this P*»3aRAH SCANS the Data pojnts of the deconvolved 
c function and determines the peak heights* location* and 
C the full width at half maximum eF EACH* AN ARTIFICIAL 
C FUNCTION li THEN BUILT FROM GAUSSIAN DISTRIBUTIONS WITH 
C THe same number of peaks* and with heights* WIDTHS* AND 
C LOCATIONS AS DETERMINED ABOVE* 

C the PEA< locations and widths are OUTPUTTED IMMEDIATELY 
C TO CHECK POR GROSS ERRORS* 

C 

dimension PEAKFN<1024I *HT(S)*LoC|S) *LR(5) 
dimension LL(5>*WF(5>*WW(9) 

READ(5D«4 iPEAKFN 
4 F0RHaT(D*1*15> 

C 

c THe following bit of code Determines the peak height* 
c location, and FWHM, 
c 

K-0 

L-0 

ISW«D 

DO 1 1-1*1024 

1F(PEA<FN( I > .LT.20* )Q0 TO I 

IF ( 9 7 4<Fn< Z ) •QT*PCAKFN( 1-1 ) *AND*PFAKFN( I ) *QT* 

#PEA<rs< Ul ) JL-L + IJ 1 Sw-i ; hT ( L ) *PEAKFN < I ) j 
4L0c(U-Ijfl0 TO 2 
1F( ISW.COtDGO TO 2 
GO TO 1 
2 <-<♦! 

IF<PEA<Fn( II .LT.HTtLI/2 * I LR ( Li-K-l j ISW-O j K-0 
1 cOMTINJC 
00 7 I-t*S 
N-LOCU I -lRI I ) 

DO I J-l*« 

1F< P£A<FN( J+N-3 I .OT.HT( I )/2. >LL< * >“L«< I 
#WF1I)..49315/MLRC1)*LL(1I)/2*>**2;Q0 TO 9 
• CONTINJC 
9 WW( 1 I -lR( IULL( I ) 
fluTPjr -Mil) 

OUTPUT LOCIII 

7 continue 

THe F3-.9WINB BIT OF CODE BUILDS THE ARTIFICIAL FUNCTION 
FROM GaUSSTanS 

DO 6 U-l*1024 

6 PEAKFNI J»-HT<ll*EXPl-lj-L0Cllll**2*WFIl) I 
A*HT(ai*EXP(-(J»L0c(2) l** 2 *taFtai I 
♦ ♦HT( 3 ) •CXPI -( J-LOC( 3 > )-*2*WF 13)) 

#*HT ( 4 ) -CXP( -I J-L0C14 ) >*»2*WFC A) I 
«-»HT( 5) -EXP(-( J-LOCIS) l»«2«wFCSI I 
WRITE! 6D*9)PEAKFN 
5 FORMAT! tOF5*0) 

STOP 
END 


li 


67 



name ■ Sj3 


C 

C this P^335*m substitutes COeFFiCIeSTS of THE AI^TIFICIaL 
c funCtI3h into the transform of t«e Deconvolved function, 

C INPUT the Humber of coefficients of the mRIGINaL 
C function that you want and the rest will Se automatically 
C substitute^ from the artificial function, 
c 

double Complex hout(io?4),mhout(io24i 
REAO( 1, 11 JHOUT 
11 FORmATj D?1 •15#2x#D2i .15) 

READj 2,22)HH0UT 
22 FORMAT) 021 .l5#2X#02l .15, 

INPUT H 

00 33 X*M+1«1025-M 
33 HOUTt I )»HHOUTt I) 

WRITE) 30.44,HflUT 

4^ FflRHAr)02l .l5.2X.D2i.i5, 

STOP 

CnO 
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NAME - P’SHA'I 


double Complex F(1024l#G(103%J«HM02«t 

READ( 8, 1 >r 

t EOR^ArO?! .19V2X«02 i*15) 

2 READI 9 , 3 )h 

3 pORiAT< 321 .1S,2X#021 .15) 

♦ Ua 5 I-1.510 

IF ( I . GT • 46 ) Q< 1 ) •Oi G < 1026- n -Oi GO TO 5 
Gt I )•( •5^4«A6,0 c0S( 3. 14 15926 53589793* 

G11025«I t«r •54+»4A*0C0S(3.i4ii92653589793* 
* I I/A5. ) ) I *IH{ 1025-1 )/F( 1025-1 ) ) 

5 CONTINJE 

NRI TE( 10. 6)G 

8 F0R*1Am2i .l5,2X#02l.l5) 

STOP 

END 


n n n n 


name . A^EA 


THIS P^33^AM FINDS THE AReAS UNOgR THE PfA<S. THE ENDPOINTS 
HUST UE OETeRHINeD BY INSPECTION. 

DIM£'^SI3N A<5)»DECON(bO30) 

DATA A/3/ 

HEA3< 73.&)DECSN 

6 F0R1AT(10F5.0) 

DO 1 I«1«a70 

1 AU ) "A( H+OECflN( 1) 

DO a I»A72#1120 

2 A(2)*A(2)4>DECnN(I) 

DO J I-l9jli2861 

3 A(3)«AO)4.0EcaN{ I) 

DO 4 I-2S16/3770 

A A( A ) >A ( A >-fDECON( I ) 

DO 5 I-AOng^AsSO 
S A ( 5 ) >A ( 5 )+OECON( I ) 

DC 7 J»l#5 

7 hRITE( 71.K)J«A(J) 

g PaR.IATC' •/ *AB£A1( *,Il, ' )-',F8.0) 

OUTPUT iThese Results must be hultiplied by the 
#PR 0PrR proportionality CHNSTaNt TO GIVE CORRECT AREA* 
STOP 
END 
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